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A DISSIPATIVE TIME REVERSAL TECHNIQUE FOR PHOTO-ACOUSTIC 

TOMOGRAPHY IN A CAVITY 

LINK V. NGUYEN AND LEONID KUNYANSKY 


Abstract. We consider the inverse source problem arising in thermo- and photo- acoustic tomography. 
It consists in reconstructing the initial pressure from the boundary measurements of the acoustic wave. 
Our goal is to extend versatile time reversal techniques to the case of perfectly reflecting boundary of 
the domain. Standard time reversal works only if the solution of the direct problem decays in time, 
which does not happen in the setup we consider. We thus propose a novel time reversal technique 
with a non-standard boundary condition. The error induced by this time reversal technique satisfies 
the wave equation with a dissipative boundary condition and, therefore, decays in time. For larger 
measurement times, this method yields a close approximation; for smaller times, the first approximation 
can be iteratively refined, resulting in a convergent Neumann series for the approximation. 


1. Introduction 

We consider the inverse source problem arising in the thermo- and photoacoustic tomography (TAT and 
PAT) |Kru95l lKru99[ [Oi^ IXW06] . It consists in reconstructing the initial pressure in the acoustic wave 
from the values of time-dependent pressure measured on a surface, completely or partially surrounding the 
object of interest |KK08| . During the last decade, significant results were obtained in solving this problem 
under the assumption that the wave propagates in free space (see, for example |FPR041 IKu07al IKuOTbl 
|Ng09llSTTO|QSUZll[IIl^lFHR071[Nan^lXWn^[PM^lP^lH£^ and reviews |KK08llKKTTllSchIl| 
for additional references). Applicability of the free space approximation depends on the type of the 
device(s) used to conduct the measurements: it is valid if reflection of waves from the detectors can be 
neglected. There are, however, a number of situations where this simple model is not applicable. For 
example, when the object is surrounded by glass plates optically scanned to measure the pressure, the 
waves experience multiple reflections as they would in a resonant cavity |CAB07[ ICB081 IKHC13| . As a 
result, the assumption indispensable in the analysis of the classical TAT and PAT about the fast decrease 
of acoustic energy within the object, can no longer be made. 

Thus, a novel approach is needed to solve the inverse source problem posed within a resonant cavity. 
The latter problem has attracted attention of analysts only recently. A particular case of a rectangular 
resonant cavity was considered in |CAB07l ICB081 IKHCISj , and several solutions that exploit symme¬ 
tries of such a geometry were proposed. In |HK15[ ISYI51IAMI5) , more general acquisition geometries 
were considered; one of the main questions investigated in these papers is the applicability of various 
modifications of the time reversal algorithm to the problem under consideration. 

Time reversal was successfully used by multiple authors, both theoretically and numerically, to solve the 
inverse source problem in the free space setting { |FPR04l IHKN081ISU091 IHri091 [QSUZl I[ ISUlll |HomI3j 1. 
In the latter case, it consists in solving the wave equation backwards in time, within the domain D 
surrounded by the detectors, using the measured data to pose the Dirichlet boundary condition. In the 
simplest case of constant speed of sound c{x) = cq, the time reversal is initialized with zero conditions 
inside the domain (at time t = T for sufficiently large T). Such approach works due to the fact that the 
solution of the direct problem in with c{x) = cq vanishes within the domain in a finite time To(f2), due 
to the Huygens principle (and so time T is chosen to be greater or equal to To(H)). In 2D and/or when c{x) 
is not constant, more sophisticated methods are used to initialize time reversal |HKN08IISU091 [QSUZl 1| ; 
however, all these techniques require a time decay of the solution of the direct problem. 

The inverse source problem in a cavity with partially reflecting walls was solved in |AM15j using 
a time reversal approach. In this case, due to the outflow of energy through the boundary of the 
domain, solution of the direct problem still exhibits time decay, and, with some modifications, time 
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reversal can still be applied. A more difficult case is that of perfectly reflecting boundary (modeled by 
zero Neumann condition on the boundary). In this case, the direct problems is energy preserving, and 
acoustic oscillations in the cavity continue forever (here we neglect the effects of attenuation of waves 
in the tissues). This immediately makes classical time reversal inapplicable: the correct pressure and 
its time derivative inside the domain are unknown at any time. The error made when replacing these 
functions by any crude guess is not small. Under the standard Dirichlet boundary condition, this error 
does not decrease as t approaches 0. 

Several attempts were made to modify the Dirichlet boundary condition so as to control this error. In 
[HK15| . the boundary condition is multiplied by a smooth cut off function rj(t/T), where rjir) vanishes at 
r = I with at least several derivatives. It was shown that under certain conditions on the eigenvalues of the 
Dirichlet and Neumann Laplacians on D, the approximate solutions corresponding to measurement time 
T converge to the correct one in the limit T —> oo. However, for an arbitrary domain the aforementioned 
conditions on the eigenvalues are difficult to verify. In |SY15) . the Dirichlet boundary condition was 
modified in such a way that the computed function is the result of averaging of a family of time-reversed 
solutions with different times T (“averaged time reversal”). The authors showed that, in the case of full 
boundary measurements, the operator that describes such a procedure is a contraction. Therefore, the 
result can be used either as a crude approximation, or a better solution can be computed by a converging 
Neumann series (involving multiple solutions of the direct problem and averaged time reversals). For the 
case when the data are available only on a part of the boundary (i.e. for the “partial data problem”), 
theoretical analysis was not completed; however, numerical experiments yielded successful reconstruction 
in this case, too. In [AM 1 5] . the problem with perfectly reflecting walls was considered mostly from a 
theoretical standpoint, and its unique solvability was proven under sufficiently general conditions. 

In the present paper we propose a new time reversal technique in which the Dirichlet boundary 
condition is replaced by a non-standard boundary condition involving a linear combination of the normal- 
and time- derivatives of the solution (see equation @). While this new condition is satisfied by the 
solution of the direct problem, the error resulting from the time reversal satisfies the wave equation with 
dissipative boundary condition, and, thus, it decreases as time t approaches 0. We call this technique a 
“dissipative time reversal”. This method can be used either directly to obtain good approximations to 
the initial pressure (which converges exponentially as T —^ oo), or as a part of a Neumann series-based 
iterative algorithmic In addition to being efficient numerically, our technique is relatively easy to analyze, 
allowing us to deliver an intuitive proof for the convergences in both full and partial data cases. 

2. Formulation of the problem and preliminaries 

2.1. The direct problem. Let D be a domain in with smooth boundary dO,. In the realistic setup of 
PAT and PAT, the dimension d equals 3. However, we will consider any d> 2, for the sake of generality. 
We will denote by u = v{x) the outward normal vector of dfl at x. The speed of sound c{x) is a positive, 
infinitely smooth function on 

For some positive time T, let us consider the following mixed boundary problem for the wave equation: 

{ Unix, f) — c^{x) Au{x, t) = 0, (x, t) € fl x (0, T], 

£u(x,t)=0, (x,t) € dn X (0,T], 

u(x, 0) = Uo(x), Ut(x, 0) = Ui(x), X € fl, 

where ug and ui are arbitrary functions with ug € Ilg(id), ui € L^(fl). When ui = 0 and uo(a:) coincides 
with the initial pressure f(x) in the tissues, the solution u(x,t) of the above problem represents a model 
of TAT/PAT within a resonant cavity formed by the reflecting walls I ICABOTT ICB081 |KHCI3j l. In this 
case, one attempts to reconstruct initial pressure f(x) from the measurements of u(x,t) made on some 
subset of boundary dfl. 

For convenience, we consider here a slightly more general problem, where both ug and ui are allowed to 
be non-zero. Let us assume that the values of pressure u(x, t) are measured on a part F of the boundary 


^ Our idea of using the Neumann series-based iterative method is inspired by the influential paper ISU09I 
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dri; we will consider both the case of full measurements when F coincides with Sil, and the case of partial 
data when F is an open proper subset of dVl. We will denote the measured data by g{z,t) : 

(2) 5 = ^^|rx[o,T]- 
Our goal is to solve the following problem: 

Problem 2.1. Find the pair {uo,ui) from g = w|rx[o,T]- dn other words, invert the map 

(3) At ■■ {uo,ui) ^ g. 

2.2. Our approach to the inverse problem. Let us fix a function A € such that A > 0 on F 

and A = 0 on 917 \ F. In addition, we will extend the data g{z, t) to dfl \ F x [0, T] by zero. The main idea 
of this paper is to find an approximation to uq and ui by solving the following time reversal problem 

( Vtt{x,t) - c^{x) Av{x,t) = 0, (a:,t) G n X [0,r], 

(4) I ■§j^v{x,t)-X{x)vt{x,t) = -X{x)gt{x,t), {x,t) e dil x [0,T], 

[ v{x,T)=0, vt{x,T)=0, X G fl. 

As we show below, for sufficiently large values of T, the functions v{x, 0) and Vt{x, 0) are good approx¬ 
imations to Uq and Ml, correspondingly. 

In particular, in the context of the inverse problem of TAT/PAT with mi = 0 and uq = /, the function 
v{x,0) is a good approximation of f{x). One can either use this approximation directly, or to realize an 
iterative refinement scheme converging to / and corresponding to a computation of a certain Neumann 
series. 

The well-posedness of problem Q, validity of a so-obtained approximation and convergence of the 
Neumann series are the subject of the following sections. The starting point, however, is the analysis of 
the direct problem 0. 


2.3. Properties of the direct problem. Let us define the space of pairs of functions H as follows 

El := {(mo,Mi)|mo G G L^(r7)}, 

where F[^(Q) is the standard Sobolev space, with the norm defined, for an arbitrary function h(x) by 


iH'i(n) 


h^{x) -f \Vh(x)\^ 


dx = 


Ii2(n) 




Then, El is a Banach space under the norm ||.|] defined by 

||(Mo,ui)f = ||MollHi(a) + l|c”^Mi||^2(n)- 

For (mo,mi) G EI, we define 

E(mo,mi) = ||VMo||i2(n) + l|c”iMi||i2(n), 

and the semi-norm |.|: 

|(mo,mi)| = [E(mo,Mi)]^/^. 

It is well known that for any (mq, mi) G EI equation Q has a unique solution u lying within the following 
class /C(r7, T): 

uGJCin,T)=C{[0,T];H\n)) n CmO,T];L^{n)). 

Therefore, At (defined by 0) is a well-defined map from EI to C{[Q,T]; 


3. Well-posedness of a mixed boundary problem 
Let A G C°°{dtd). In this section, we study the following mixed boundary value problem: 
wtt{x, t) — c^{x) Aw{x, t) — 0, (x, t) gVLx (0, T), 

^w{x, t) + X{x) wt{x, t) = htix, t), (x, t) G drix (0, T), 

w{x,0) = wo(x), Wt{x,0) = wi{x), xGVL. 


( 5 ) 




4 


LINK V. NGUYEN AND LEONID KUNYANSKY 


Since this is not a standard problem, let us first define its weak solution. We will follow the spirit of 
[BLR92| . where the weak (distributional) solution is defined for the case h = 0. For our purposes, it is 
sufficient to assume {wo,wi) € El and h S C([0, T]; 

Definition 3.1. Given {wo,wi) S El and h G C{[0,T]; (dfl)), we say that w G x (0,T)) is a 
solution of if the following equation holds for any test function if G x (0, T)): 

{Wjtjj) = — // [z(;o(a;)0) — 0)] dx + / X{x) Wq{x) 'i>{x,0) dx 

JJnx[0,T] JdQ 

(6) ~ h{x,t)'^t{x,t) dx dt — / /i(a;, 0) ^'(a;, 0) dx. 

d7anx[o.r] Jan 

Here, is the solution of the dual problem 

( c~'^{x)'i>tt{x,t) - A'i>{x,t) = ip{x,t), {x,t) G XI X {0,T), 

< d,^'i>{x,t) — X{x)'i’t{x,t) = 0, {x,t) G on X {0,T), 

[ 4'(x, T) = 0, T) = 0, xGQ. 

The uniqueness of solution w is obvious, since {w, ip) is defined by the right hand side of Q for all test 
functions ip. The existence and regularity of this solution are more complicated. We only present here a 
partial result that is needed in the present paper: 

Proposition 3.2. Assume that h = 0 and {wo,wi) G H. Then, problem has a unique solution 

w G 1C{n,T). 


Let us recall a result by irEro] . Assume that h = 0 and {wo,wi) G id^(Sl) x H^(fl) satisfies the 
compatibility condition 

d 

(7) ~ 

Then, problem ([^ has a unique solution 

w G C([0, T];H^(fl)) n Ci([0, T];H^(n)) n C^([0, T];L^(fl)). 


Proposition 3.2 can be proved by a simple approximation argument which we present here for the sake 
of completeness. 


Proof of Proposition \3.2\ It suffices to prove the existence part, since the uniqueness is trivial as 
mentioned above. Consider {ipQ,(pi) G H^{Xl) x H^{Xl) such that: 

1) iTo,Ti) -t (wqjWi) in H^{n) x L'^{Xl), 

2) d,,ip^\dn = 0 and 'Pi\an = 0. 

Such a sequence {{^Pg, Ti)} always exists. For instance, we can choose to be a linear combination 
of eigenvectors of the Neumann Laplacian. Meanwhile, pAp can be chosen as a linear combination of 
eigenvectors of the Dirichlet Laplacian. 

Consider problem (§ with the initial condition (instead of {wq,Wi)). Since ((^q,:/?") satisfies 

the compatibility condition ([^, problem ([^ with such initial condition has a unique solution 

w” G c([o, T\,H‘^iyt)) n c^([o, t]-,h^{xi)) n c‘^{[q,T]-,L^ iyt)). 

Moreover, by simple integration by parts, we obtain for any fg G 


Therefore, 


E(w”-u;'",to)+ f A(x)K-u>)"pdxdt = E(u>”-u;™,0). 

J 30x [0,io] 

E(zc" - < E(u;" - u;'",0). 


E(u;" - u;-,0) = ||V(^^ - + Wc-^T^i - tT 


IlRO) 


0 . 


We notice that 
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E(k;" - w^,t) = II .) - .))f + ||c-i«(i,.) - .))f, 

we obtain {w"} is a Cauchy sequence inC([0, T];H^{Q)) and {w^} is a Cauchy sequence inC([0, T]; L^(f2)). 
Therefore, there exists a function w G IC{il.,T) such that 

{w”} -^win C{[ 0 ,T],H\n)), and {<} ^ Wt in C{[ 0 ,T]; L^{n)). 


It is easy to verify that rc is a 
the existence. 


solution of problem ([^, from the definition 3.1 


This finishes our proof for 

□ 


Let us apply the above result to the time reversal problem @- 
Proposition 3.3. Problem 0 has a unique solution v G /C(^2,T). 

Proof. We first notice that problem Q is the time reversed version of problem ([^, considered in Section]^ 
The uniqueness of the solution follows trivially from the definition |3.1| We now prove the existence. Let 
us consider the following problem 


r Uttix,t) - c^{x) AU{x,t) =0, (x.t) G n X [0,T], 

(8) i ^U{x,t) - X{x)Ut{x,t) = 0, (a;,f) G X [0,r], 

[ U{x,T) = u{x,T), Ut{x,T) = ut{x,T), x G it. 


where u{x,t) is the solution of the direct problem Q. Since {u{.,T),ut{.,T)) G H, we obtain from 
Proposition 3.2 that the above problem has a (unique) solution U G 
We notice that u also satisfies the boundary condition in problem @- 

Let V = u — U. It is easy to verify that is a solution of Moreover, from the regularity of u and 
U, we conclude that v G /C(n, T). □ 


4. Solution of the inverse problem 


4.1. Contraction properties of the time reversal operator. Similarly to |SY15lfSM15] . our analysis 
on the inversion of A relies on known results on stabilization of waves |BLR92| . For the sake of simplicity, 
we assume that all the geodesics of (M^, c~^ dx^) have finite contact order with the boundary Under 
this condition, the generalized bi-characteristics of the wave operator U = dtt — (a:) A on U are uniquely 

defined (see, e.g., [BLR,92| L Their projections on the physical space (i.e., U) are called the generalized 
rays. 

Throughout the paper, we will assume that the following condition is satisfied: 

Condition 4.1. There is a finite value T{n,T) > 0 such that every generalized ray of length^T{^l,T) 
intersects T at one (or more) non-diffractive point(s). 


This is the geometric control condition (GCC), well known in control theory. We refer the reader 
to [BLR92] for a detailed discussion of GCC. This condition was shown in |SY151 IAM15] to be sufficient 
for the stability of the inversion of A. Our results are also based on the assumption that GCC is satisfied. 


Our work is based on the following fundamental result due to [BLR92] : 

Proposition 4.2. Consider problem 0 with (7 = 0. Assume that Condition \).l\ holds. Then, there is 

\iui.,T),Uti.,T))\ < S{T)\{uo,Ui)\. 


S{T) < 1 such that 


Remark 4.3. By applying Proposition \4.^ k times, one obtains 

\{u{.,kT),ut{.,kT))\ < S'^{T)\{uo,u,)\. 

It is easy to conclude that \{u(.,T),ut{.,T))\ -G 0 exponentially as T ^ oo. That is, there are constant 
(71(12) and a > 0 such that 

\{ui.,T),Ut{.,T))\ < Ci(U)e-“^ |(uo,ni)|. 


^ The length here is understoo(d in the metric c ^ dx^. Condition |4.l| means that all the singularities of the solution of 
the wave equation 0 that start propagating at time 0, traveling along generalized bi-characteristics, reach the set T within 
the time interval [0,T(Q,r)]. 
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Given boundary data g{x, t), we define the time reversal operator A through the solution v{x, t) of the 
problem Q: 

Ag = (n(.,0),ut(.,0)). 

Lemma 4.4. Assume that condition ED holds and T > r(n,r). Let I denote the identity map, then 
there is S{T) < 1 such that 

\{I - AAt){uo,ui)\ < (5(r)|(uo,Mi)|, for all {uo,ui) G H. 

In other words, I — AAt is a contraction on IHI under the semi norm | • |. 

Proof. Let (mo,mi) G EI and u be the solution of Q. We observe that 

(/-AAt)(uo,ui) = ([/(., 0), [/*{., 0)), 

where U is the solution of the problem ([^. 

Due to the well-known conservation of energy in the solution of the wave equation with Neumann 
boundary condition, 

E(uo,ui) = E(M(.,T),ut(.,T)), 
or 

|( mo , mi )| = \iu(.,T),uti.,T))\. 

Due to Proposition |4.2[ 

\{U{.,0),Ut{.,0))\ < S{T)\{u{.,T),ut{.,T))\, 

and the proof follows. □ 


Remark 4.5. Using Remark \4-3[ instead of Proposition 14-^ in the above proof, we obtain that there 
exist constants (71(0) and a > 0 such that 

|(/ - AAkT)iuo,ui)\ < ( 7 i( 0 )e““^|(uo, 1x1)1. 

That is, the induced seminorm of I — AAj- decreases exponentially as T -^ 00 . 


Lemma 4.4 describes the contraction property of our time reversal operator in the semi-norm |.|. By 
itself, this result is not sufficient to prove the convergence under the norm ||.||. However, by restricting 
our attention to appropriate subspaces of H, such a convergence can be attained. Indeed, let us consider 
the subspace Ho of H defined by 


= < h = {ho hi) G 


hodx = 0 > , 


an 


and the subspace Hi of Hq consisting of pairs with the second component equal to zero: 

Hi = {h = (hoO) G Hq} . 

Let us introduce two projectors Hq and Hi mapping the elements of H into Hq and Hi, respectively: 


Hoh = Ilo{hohi) = I /iQ — J ho, hi j , 

V an J 

Hih = ni(/io,/ii) = I ho — J ho, 0 I , 

\ an J 


where |cID| is the surface area of dUl. These projectors do not increase the semi-norm |.|. That is, 

|noh|<|h|, |nih|<|h|. 

Moreover, the subspaces Ho and Hi are invariant under compositions Hg)/ — AAt) and Hi)/ — AAt). In 
addition, 


Vh gHo, 
Vh G Hi. 


no(/ — AAt)\x={I — HoAAt")!!, 
ni(/ - AAT)h={I - HiAArjh, 
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Therefore, in accordance with Lemma 4.4 operators (/ — IIoAAj’) and (/ — IfiAAr) are contractions in 
Hq and Hi correspondingly, under the seminorm |.| : 

(9) |(/-noAAT)h|< |(/-AAT)h| <5(T)|h|, Vh G Hq, 

(10) |(/-niAAT)h|< 1(7-AAT)h| < ,5(T)|h|, Vh G Hi. 

Now we can invert operator Ay, restricted to Hq, by constructing converging Neumann series 

OO OO 

^(/ — IIoAAt)^IIqAAt{uo, ui) — — IloAAT)^IloAg. 


/c=0 


k—O 


Below we show that this series converges not only under the semi-norm |.| but also in the norm ||.||. In 
other words, we claim that the partial sums u^”)of these series 

k 

(11) u('=) = ^(7 - n^AkTyn^Ag 

j=o 

converge to u = (mo,mi) in the norm ||.||. These partial sums can be easily computed by the following 
iterative algorithm 

u(o) = 0, 

(12) = {I- noAAT)u('=) -1- UoAg. 


The following theorem gives us a solution to Problem |2.1| by a Neumann series: 

Theorem 4.6. Suppose that condition \4-l\ holds and the observation time T satisfies T > T(n,r). Then, 
the iterations defined by (11) (or, equivalently, by ll^) converge to u in norm ||.||, as follows: 

|lu-u('=)|| < Cp(fl)5(T)'=||u||, for all ugUo, 

with some constant Cp{VL) > 1 specified below. 

To prove the above theorem, we will need the following generalization of the Poincare inequality. 

Lemma 4.7. There is a constant Cp{n) > 1 such that for all h G 77^(17) satisfying Jg^hdx = 0 the 
following inequality holds: 

Proof. Indeed, assume that the above statement is not true. Then, there exists a sequence {/in})0Li C 
77^(17) such that 

l|/in||L 2 (o) = I and lim ||V/i„||p 2 (n) = 0. 

It follows that {hn} is bounded in 77^(17). Therefore, there is a weakly converging in 77^(17) subsequence 
{^nk}T=i function u G 77^(17), such that, 

hn,, —t u strongly in 

huk u weakly in 77^(17), 

hnk —t u weakly in 

The limit u has the following three properties 

lkllLVO) = l, l|Vu||L 2 (f 2 ) = 0, udx = 0. 

an 

The last two equations yield u = 0, which contradicts to the first property ||m||l 2 (q) = I, thus completing 
the proof. □ 

Corollary 4.8. There is a constant Cp(fil) > 1 (given by the above Lemma) such that the following 
inequality holds: 


(13) 


|h| < ||h|l < C'p(f7)|h|, forallhGUo. 
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Proof of Theorem \4-€\ By a simple induction argument applied to the recurrence relation ( |12[ ) one 
obtains the following identity: 

u - = (/ - noAAT)'=u, A: = 0,1, 2,3,... 

Since u G Ho, applying inequality (|^ results in the following inequality 

<<5(r)'=|u|. 


Further, using inequalities (13), one can transition to the following norm estimate 

< Cp{n)5{Tf' 


lu-u('=)| 


U 


□ 


thus, proving convergence of the Neumann series in the norm ||.||. 

Remark 4.9. We note that, due to Remark \4-5\ 

(14) Iju-noAffll <Ci(F!)Cp(fl)e-“^||u||. 

Therefore, when T is sufficiently large, = PoAg is a good approximation to u. 

4.2. Inverse problem of TAT/PAT. The inverse problem of TAT/PAT is a particular case of Prob¬ 
lem with Uq = f and Ml = 0. Therefore, / can be recovered from g by using the algorithm described 
in the previous section. However, the convergence can be accelerated and computations simplihed by 
projecting computed approximations onto space Hi rather than Hq. 

Let us introduce the notation f = (/, 0). The measured data g are still defined by equation ([^ with 
u{x, t) being a solution of the direct problem ([^ with initial conditions 

(M(0,a:),Mt(0,a:)) = f(a;). 

We show below that the partial sums 

k 

uW = ^(/ _ liiAKTyWAg- 
1=0 

of the following Neumann series converge to f in norm ||.||. These sums are easy to compute using the 
following iterative relation 

u(o) = 0, 

(15) = {I- niAAp)u('=) AiAg. 


Theorem 4.1 0. Assume that the initial conditions of Problem 2.1 is given by (mq, Mi) = f. Suppose also 
condition 4-1 *5 satisfied and T > T(n,r). Then, iterations defined by (11) (or, equivalently, by 


(12)) converge to f in norm ||.||, as follows 

llf-u«||<Cp(fl)5(T)'=||f||. 

Proof. The proof is almost identical to that of Theorem 4.6 with inequality (10) used instead of ([^. □ 

Remark 4.11. Similarly to Remark \4.9l 

\\i -U^AgW < cyn)Cp{n)e-^^\\f\\, 

and, when T is sufficiently large, = AiAg is a good approximation to f. 

5. Numerical implementation and simulations 

5.1. Implementation. One of the advantages of the present method is the ease of implementation 
using standard finite differences. Unlike algorithms of [SY15] , our approach does not require solving 
the Dirichlet problem for Laplace equation to initialize the time reversal. (The latter problem is well 
studied and various methods for its solution are known. However, efficient numerical schemes for arbitrary 
domains are quite sophisticated and require noticeable effort to implement). 

Our numerical realization of the algorithm is based on equation (15) that requires computing operators 
Ap and A. These operators represent solutions of the wave equation forward and backwards in time, 
respectively; they were calculated using finite difference stencils as described below. 
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(a) Phantom (b) Full data reconstruction (c) Partial data reconstruction 



Figure 1. Reconstruction with T = 5. In (d): gray line represents the phantom, dashed 
line shows image reconstructed from the partial data, black line shows full data recon¬ 
struction 


Our simulations were performed on a 2D square domain [—1,1] x [—1,1]. Throughout this section we 
will represent our 2D spatial variable x in the coordinate form, and will change the notations for all 
functions correspondingly: 


X = (x, y), u{x, t) = u(x, y, t), v{x, t) = p(x, y, t), etc. 


Our square domain was discretized using Cartesian grid of size 257 x 257 , with the step Ax = Ay = 2/257; 
time was discretized uniformly with the step At = 0.5Ax. The speed of sound c(x,y) was set to 1, for 
simplicity. Time stepping inside the domain in the forward direction was implemented by applying 
standard second-order centered stencils in both time and space to the discretized solution u{xk,yi,tj) 


02 


92 




9y2 


92 


u{xk,yutj) w ^u{xk,yi,tj) = 


L{xk+i,yi,tj) + u{xk-i,yi,tj) - 2u{xk,yutj) 

Ax2 

_ u(xfc,y/+i,tj) -I- u{xk,yi-i,tj) - 2u{xk,yi,tj) 


Ax2 


r u{xk,yi,tj+i)+ u{xk,yi,tj_i)-2u{xk,yi,tj) 

^u{xk,yi,tj) w ^u{xk,yut,) = --, 


resulting in the formula 


92 


92 


M(xfc,y/,tj+i) = 2u{xk,yi,tj) - u{xk,yi,tj-i) + Ar ( —u{xk,yutj) + —u{xk,yi,tj) ) , 


9y2 
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(a) Exact data g{x,t) and data with added 50% noise (in norm), for one x 



(b) Phantom (c) Full data reconstruction (d) Partial data reconstruction 



Figure 2. Reconstruction with T = 5 and with added 50% noise (in norm). In 
(e): gray line represents the phantom, dashed line shows image reconstructed from the 
partial data, black line shows full data reconstruction 


where tilde denotes approximate quantities. Time stepping backwards in time (when computing A) were 
done similarly 


/ Q2 Q2 

(16) - vixk,yiAj+i) + ( -^v{xk,yi,tj) + 

When computing action of the operator At (forward problem), the Neumann boundary condition was 
represented by the simplest first-order two-point stencil; this results in the values at the boundary points 
being set to the values of the nearest grid points. 
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The discretization of the non-standard boundary condition 
(17) 


d d d 

— v{x,t) - A(x) -^^v{x,t) = -\{x) —g{x,t), 


arising in problem Q was performed as follows. The simplest hrst-order two-point forward stencils were 
used to approximate all the derivatives; for example was approximated at t = tj_i, x = xq by 

d , . \ 9 f _ y{^o,yhtj) - v{xo,Yi,tj_i) 

and was computed similarly. The normal derivative was also approximated by the simplest two-point 
stencil applied to values at time tj_i ; for example, one the side with x = xq the following formula was 
used: 

—v{xQ,Yi.,tj-i) = -—v{xo,Yhtj-i) w -—v{xo,Yi,tj-i) = --. 

Substituting the last two equations into 0 resulted in 

v{xQ,Yi,tj-i) - v{xi,Yi,tj-i) _ X{xq,yi) 


Ax 


At 


b(xo, Yhtj) - w(xo, Yhtj-i) - 5(xo, Yhtj) + ff(xo, 


or 

(18) 


K^o,Yi,tj-i) = ^'(xl,y/,^J-l) +7(xo,y0 b(xo,yz,ti) - v{xo,Yhtj-i) -5(xo,yz,ii) +5(xo,y/,tj-i)] 


where 


A(xo,yz)Ax 
7(xo,yz)^--■ 


Solving (18) for ?;(xo,y/, tj_i) yielded 


(19) 


^'(xo,y/,tJ-l) = 


v(xi,Yi,tj-i) 7(xo,y;) [vi^o,Yhtj) -5(xo,y/,ij) + 9i^0,Yi,tj-i)] 


l + 7(xo,y;) l + 7(xo,yi) 

Approximation of boundary condition 0 on other parts of the boundary was done similarly. In order 
to apply (19) (and similar expression on the other parts of the boundary), one first applies (16) at all 
discretization points inside of the computational domain. Then (19) is fully defined. 

In the absence of experimentally measured data, in order to validate the reconstruction algorithm one 
needs to simulate values of g{x,t) on T. One could use the finite difference algorithm described above to 
approximately compute g{x,t) for a chosen phantom f = (/, 0). However, doing so would constitute the 
so-called “inverse crime”: sometimes simulations will produce inordinately good reconstructions due to the 
spurious cancellation of errors if the forward and direct problems are solved using the same discretization 
techniques. Thus, in order to compute g we utilized the following method based on separation of variables. 
Function / was expanded in the orthogonal series of eigenfunctions (fk,i of the Neumann Laplacian on 
our square domain 

/(x,y) = ^Cfc,zV5fc,z(x,y), 

k,l 

<Pfe,;(x,y) = cos(fcx) cos(Zy), /c = 0 , 1 , 2 ,..., ; = 0 , 1 , 2 ,..., 

X = 7r(x -b l)/2, y = 7r(y -b l)/2. 


This was done efficiently using the 2D Fast Cosine Fourier transform algorithm (FCT). Then, solution of 
the forward problem was computed as the series 

?/(x,y,t) = ^Cfc,/i^fc,/(x,y)cos(Afc,/i), Xkj = ^ kj = 0, 1 , 2 ,.... 

k,l 

For each value of t, the above 2D cosine series were also summed using the FCT. The resulting algorithm 
is quite fast, and, more importantly, it is spectrally accurate with respect to /. If / has high degree 
of smoothness, this series solution yields much higher accuracy than the finite difference techniques we 
utilized as parts of the reconstruction algorithm. 
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(a) Phantom (b) Initial approximation IViAg 



(c) Iteration #2, (d) Iteration ^5, 



Figure 3. Reconstruction from the full boundary data, T = 1.6. In (e): gray line is the 
phantom, dashed line shows the initial approximation, black line presents iteration ^5 


5.2. Simulations. We conducted several numerical experiments to verify theoretical conclusions of pre¬ 
vious sections. Two acquisition schemes were considered: a full data scheme with g{x, t) given on all sides 
of the square domain, and a partial data scheme with g known only on the left and bottom sides of the 
square. Since we assumed c{x) = 1, T{fl, P) equals 2^/2 in the former case (i.e., the length of the diagonal 
of the square) and 4\/2 in the latter case. Experimentally, we found that these times are too pessimistic, 
and a half of that time is quite enough for the reconstruction. This implies that our theoretical results 
may be not sharp. Although we were not able to improve these theoretical estimates, we present below 
simulations with the measurement times significantly reduced compared to T(fl,r). 
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(a) Phantom (b) Initial approximation IViAg 




(e) Central horizontal cross-sections of (a), (b), and (d) 


Figure 4. Reconstruction from the data given at the left and bottom sides of the 
square domain, T = 3. In (e): gray line is the phantom, dashed line represents the initial 
approximation, black line shows iteration ^5 


As a phantom, we utilized a sum of six shifted finitely supported (7^(17) functions of the radial variable, 
shown as a color-scale image in Figure 1(a). Such a smooth phantom reduces errors related to finite 
difference computations and allows us to concentrate on convergence of the algorithm per se. 

The goal of our first two simulations was to see how well the initial approximation TiiAg can approx¬ 
imate /. To this end the measurement time T was chosen to equal to 5; this corresponds to the time 
of two and a half bounces of waves between the opposite sides of the domain. Reconstruction from 
full data is shown in Figure 1(b) and partial data reconstruction is presented in Figure 1(c). Image in 
Figure 1(d) demonstrates profiles of the central horizontal cross sections of the three previous images. 
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with the phantom represented by a gray line, full time reconstruction shown as a black line, and partial 
data image drawn with a dashed line. The full data image is practically perfect; the corresponding black 
line in Figure 1(d) lying on top of the gray line, rendering the latter almost invisible. The relative 
reconstruction error equals 1.1% in the image of Figure 1(b). The partial data reconstruction is just 
slightly less accurate, with the relative Lp‘{Vl) error equal to 6.8%. In any reconstruction from real data 
such error would be negligible compared to errors introduced by imperfections of real data. 

In order to illustrate the low sensitivity of the algorithm to the noise in the data, we repeated the 
previous simulation with the data contaminated by 50% white noise (in the relative LF' norm). As in the 
first simulation, only the initial guess I^iAg was computed, without the successive iterative refinement. 
Figure 2(a) shows the time series representing g{x,t) for one of the points x, with and without added 
noise. Figure 2(b)-(d) demonstrate the same phantom, and the full and partial data reconstructions, 
respectively. The errors in the two reconstructions were of about the same order, 19% and 22% in the 
relative norm, with the partial data giving, for some reason, slightly better result in this norm. 

The central cross sections of these images are shown in Figure 2(e). 

Our remaining two simulations were intended to verify the convergence of the algorithm in the case 
when measurement time T is close to a half of T(il,F). Figures 3(a)-(e) demonstrate results of a full 
data reconstruction with T equal to 1.6 (compare to T{V,,T) = 2\/2 w 2.828). Figures 3(a)-(d) show 
the phantom, the first approximation IIiAg, and the second and the fifth iterations and u^®^), 

respectively. Figure 3(e) presents central horizontal cross-sections of the phantom, the first approximation 
WiAg, and of the fifth iterations. One can notice that, while the initial approximation had been noticeably 
distorted, the fifth iteration yields a close approximation to f{x). The relative L^(fl) norm of the error 
in was 3.3%. 

The final series of images in Figure 4 demonstrates the results of the reconstruction from the partial 
data, with T equal to 3 (compare to T(n, F) = 4-\/2 w 5.6569). As before, the phantom, the initial approx¬ 
imation, and the second and the fifth iterations are shown in Figure 4(a)-(d), respectively. Figure 4(e) 
presents the central horizontal cross sections of images in Figure 4(a), (b), and (d). The relative L^(n) 
error in the fifth iteration was 5.4%; for most practical purposes this would be more than acceptable. 


6. Conclusions 


We presented a novel dissipative time reversal approach for solving the inverse source problem of 
TAT/PAT posed within a cavity with perfectly reflecting walls. Unlike the previous work of |HK15l[SY15) 
where Dirichlet boundary condition was used, we utilize the non-standard boundary condition Q. The 
latter leads to the dissipative boundary condition (§ imposed on the error U{x,t), and, hence, to a 
natural decay of U{x, 0) with the growth of T. Our approach results in two reconstruction methods: i) a 
non-iterative approximation, converging exponentially to / as T —>■ oo, and ii) a Neumann series formula. 
These two algorithms are applicable for both full and partial data problems. 

Compared to |HK15] , where rather stringent conditions on the eigenvalues of the Neumann and Dirich¬ 
let Laplacians on U are required for convergence, our approach is based on the much less restrictive GCC 
(Condition |4.1[ ). Moreover, unlike the method of |SY15| . our technique does not require computing the 
harmonic extension of the boundary values, which significantly simplifies its numerical realization. It 


should be noted that the requirement T > r(U,F) is sharp for the convergence in Theorem 4.6 (dealing 
with the general Problem |2.1[ ). However, it is not sharp for the convergence in Theorem 4.10 (dealing with 
the inverse problem of TAT/PAT). Specifically, for the case of the full data, it is twice of the sharp time 
for the convergence of the TAT/PAT’s inverse problem, obtained in [SY15| . Nevertheless, our numerical 
simulations show that the present algorithm performs very well with such sharp measurement times. 

While we only considered the simplest wave equation in the Euclidean spaces, our analysis can be 
extended to problems formulated on Riemannian manifolds (as in [SY15| 1 and/or to problems with a 
potential (as in [AM15] b 
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